Data

First we import our tidied salamander morphology & trait dataset.

# Import tidied salamander eye size and trait data
salamanders <- data.frame(read.csv("../Data/Tidy/salamanders_tidy.csv", header=TRUE, na.strings=c("", "NA", " ","<1")))

# Quick look at data structure
str(salamanders)

Next, we import the amphibian tree from Jetz and Pyron (2019) that has been pruned to match the salamander species in this dataset.

# Import pruned phylogeny
caudatatree <- read.nexus("../Data/Tidy/caudata-tree.nex")

# Plot tree
plot.phylo(caudatatree, show.tip.label=TRUE, cex = 0.7)

Data exploration

First, we can take a preliminary look at all of our data so that we can see how it looks and whether we notice any obvious outliers that we might want to go back and check out. Here we plot out our measurements from the right side vs. the left side of each specimen so that we can make sure that we have consistent measurements throughout the dataset. We use the package plotly to create interactive plots - hover over a point to see the name of the species measured.

Right eyes vs. left eyes

We expect that measurements of left and right eyes and corneas should be similar within specimens. We can plot out our data to test for this and make sure we don’t see any glaring issues. We can also fit a regression, which we would expect to have a slope near 1 and an intercept near 0 if left and right eye measurements are similar.

#### right eye diameter vs. left eye diameter ####

# Linear model 
RLeye.lm <- lm(ED_right_mm ~ ED_left_mm, data = salamanders)
summary(RLeye.lm)
## 
## Call:
## lm(formula = ED_right_mm ~ ED_left_mm, data = salamanders)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.56832 -0.13433 -0.01103  0.13191  0.67936 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.016679   0.035501    0.47    0.639    
## ED_left_mm  1.009053   0.009215  109.51   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2063 on 391 degrees of freedom
##   (45 observations deleted due to missingness)
## Multiple R-squared:  0.9684, Adjusted R-squared:  0.9683 
## F-statistic: 1.199e+04 on 1 and 391 DF,  p-value: < 2.2e-16
# Plot
RLeye.plot <- ggplot(salamanders, aes(x = ED_left_mm, y = ED_right_mm, text = Genus_Species)) +
 geom_point(alpha = 0.9) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  geom_abline(slope = coef(RLeye.lm)[[2]], intercept = coef(RLeye.lm)[[1]]) 


# Interactive plot
ggplotly(RLeye.plot)

As we would expect, measurements of the left and right eye are highly correlated, with a slope close to 1 and an intercept close to 0. We can do the same test for the corneal measurements.

#### right cornea diameter vs. left cornea diameter ####

# Linear model 
RLcornea.lm <- lm(CD_right_mm ~ CD_left_mm, data = salamanders)
summary(RLcornea.lm)
## 
## Call:
## lm(formula = CD_right_mm ~ CD_left_mm, data = salamanders)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.70444 -0.11965 -0.01275  0.14285  0.79425 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.11282    0.03227   3.497  0.00052 ***
## CD_left_mm   0.97073    0.01163  83.458  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.2149 on 434 degrees of freedom
##   (2 observations deleted due to missingness)
## Multiple R-squared:  0.9413, Adjusted R-squared:  0.9412 
## F-statistic:  6965 on 1 and 434 DF,  p-value: < 2.2e-16
# Plot
RLcor.plot <- ggplot(salamanders, aes(x = CD_left_mm, y = CD_right_mm, text = Genus_Species)) +
 geom_point(alpha = 0.9) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  geom_abline(slope = coef(RLcornea.lm)[[2]], intercept = coef(RLcornea.lm)[[1]]) 


# Interactive plot
ggplotly(RLcor.plot)

Again, we find that left and right corneal measurements are highly correlated, with a slope close to 1 and an intercept close to 0. Together, these results indicate that measurements were consistent between the left and right eyes of a specimen, which is good! For looking at interspecific allometry, we will calculate species means from these data so that we can use phylogenetic comparative methods.

Interspecific allometry

Here, we will use phylogenetic generalised least-squares (PGLS) regressions to see how eyes are scaling with body size across different species of salamanders with different ecologies. PGLS analyses are similar to normal linear regressions, but also account for the non-independence of data collected from close relatives due to shared evolutionary history.

In these analyses, we will log eye size and body size variables, as this is standard for allometric relationships. First, we will calculate the mean measurements across specimens for each species. We need to make sure in doing this that we omit rows with missing data for a variable of interest. In looking at the data, we can see that we have many more measurements for cornea diameter than for eye diameter.

Species means for morphological comparisons

We then find species means for various morphological parameters. There are some specimens that lack data for eye size but have data for svl, mass, and corneal diameter (these measurements are present for all specimens). However, when eye size measurements are missing, they are consistently missing for all specimens measured within a species. This means we don’t have to worry about missing data for one specimen corrupting the mean across several specimens for eye size (it is either present for all specimens in a species or absent for all specimens in a species). So we can make just one dataframe of species averages, and some species will just lack eye size data.

# Species means for eye size, cornea size, SVL, and rootmass data
av.sal <- salamanders %>% 
  mutate_if(is.character, as.factor) %>%
  group_by(Genus_Species) %>%
  summarise(eye_av = mean(eyemean), 
            cor_av = mean(cormean),
            svl_av = mean(SVL_mm), 
            rootmass_av = mean(rootmass), 
            n = n())

# Merge data means with other column info, keeping only trait columns
av.sal <- merge(av.sal, salamanders[match(unique(salamanders$Genus_Species), salamanders$Genus_Species), ], by="Genus_Species", all.x = TRUE, all.y = FALSE) %>% 
  select(Genus_Species, eye_av, cor_av, svl_av, rootmass_av, n, Order, Suborder, Family, Genus, Species, Gill_Presence, Activity_Period, Adult_Habitat, Development, Life_History, Larval_Habitat) 

# Remove Eurycea rathbuni and Proteus anguinus (blind salamanders) from dataframe
av.sal <- av.sal[-c(77, 133),]
# Export tidied dataset for analyses
write.csv(av.sal, file = "../Data/Tidy/salamanders_averages.csv")

Color figures with pallettes in the frog paper

Here we run a custom color palette for the adult habitat categories to be used for all following figures.

First we define what we want that color palette to be:

# See what the levels of adult habitat are named
levels(as.factor(av.sal$Adult_Habitat))
## [1] "aquatic"      "scansorial"   "semiaquatic"  "subfossorial" "terrestrial"
# Define a colorblind-friendly vector of colors for adult habitat
col_hab <- c("aquatic" = "#0072B2",
             "scansorial" = "#009E73",
             "semiaquatic" = "#56B4E9",
             "subfossorial" = "#CC79A7",
             "terrestrial" = "#E69F00")

# Now you can see that each habitat is assigned a hex color when you look at the vector
col_hab
##      aquatic   scansorial  semiaquatic subfossorial  terrestrial 
##    "#0072B2"    "#009E73"    "#56B4E9"    "#CC79A7"    "#E69F00"
# Make a vector of point shapes for adult habitat
shape_hab <- c("aquatic" = 25,
             "scansorial" = 24,
             "semiaquatic" = 15,
             "subfossorial" = 16,
             "terrestrial" = 18)
#see that each state has a shape now
shape_hab
##      aquatic   scansorial  semiaquatic subfossorial  terrestrial 
##           25           24           15           16           18

Phylogenetic generalised least squares regressions

In order to fit a PGLS regression in caper, we first need to make a comparative data object that includes our dataset and our phylogeny. Note that the gilled Pleurodeles waltl will be dropped from this, as there can only be one row for each species.

# Make row names of dataset the species names (so it will match phylogeny tips)
rownames(av.sal) <- av.sal$Genus_Species

# Check that names match in dataframe and tree
name.check(phy = caudatatree, data = av.sal)

# Use caper function to combine phylogeny and data into one object (this function also matches species names in tree and dataset)
sal.comp <- comparative.data(phy = caudatatree, data = av.sal, 
                            names.col = Genus_Species, vcv = TRUE, 
                            na.omit = FALSE, warn.dropped = TRUE)

# Check for dropped tips or dropped species
sal.comp$dropped$tips #phylogeny
sal.comp$dropped$unmatched.rows #dataset

Here we fit PGLS regressions and plot the resulting outputs for the logged eye size and body size variables.

PGLS of log eye diameter vs. log cube root of mass

# PGLS eye diameter vs. cube root of mass
pgls_edmass <- pgls(log10(eye_av) ~ log10(rootmass_av), 
               data = sal.comp, 
               lambda = "ML", #uses Maximum Likelihood estimate of lambda
               param.CI = 0.95)

# Diagnostic plots
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_edmass)

par(mfrow = c(1, 1))

# Likelihood plot for Pagel's lambda. Solid red line indicates estimate for lambda and broken red lines indicate the 95% confidence interval
lambda.edmass <- pgls.profile(pgls_edmass, "lambda")
plot(lambda.edmass)

# Print model output 
summary(pgls_edmass)
## 
## Call:
## pgls(formula = log10(eye_av) ~ log10(rootmass_av), data = sal.comp, 
##     lambda = "ML", param.CI = 0.95)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0283788 -0.0042682  0.0000067  0.0040428  0.0180521 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.845
##    lower bound : 0.000, p = 6.8369e-11
##    upper bound : 1.000, p = 9.5763e-07
##    95.0% CI   : (0.654, 0.943)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                    Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)        0.392719   0.042518  9.2365 4.441e-16 ***
## log10(rootmass_av) 0.798584   0.043509 18.3546 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.00678 on 137 degrees of freedom
##   (16 observations deleted due to missingness)
## Multiple R-squared: 0.7109,  Adjusted R-squared: 0.7088 
## F-statistic: 336.9 on 1 and 137 DF,  p-value: < 2.2e-16
# Quick plot of pgls
plot(log10(eye_av) ~ log10(rootmass_av), data = av.sal)
abline(pgls_edmass)

Plot of eye diameter versus cube root of mass

# Plot average eye diameter (mm) against cube root of mass (g)
plot_pgls_edmass <- ggplot(av.sal, aes(y=eye_av, x=rootmass_av, text=Genus_Species)) +
  geom_point(alpha = 0.9, aes(color = Adult_Habitat), size = 5, na.rm = TRUE) + 
  scale_color_manual(values = col_hab, name = "Adult habitat") +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(size = 26),  
    axis.text.y = element_text(size = 26), 
    axis.title.x = element_text(size = 30, color = "black", margin = margin(t = 10)), 
    axis.title.y = element_text(size = 30, color = "black", margin = margin(r = 18)), 
    axis.ticks.length = unit(0.3, "cm"), 
    panel.border = element_rect(size = 2)  
  ) +
  scale_y_log10(name = "Eye diameter (mm)", breaks = c(2, 3, 7)) +
  scale_x_log10(name = "Cube root of mass (g)", breaks = c(1, 3, 5), limits = c(0.5, 5.5)) +
  geom_abline(slope = coef(pgls_edmass)[[2]], intercept = coef(pgls_edmass)[[1]], linetype = "solid", size = 1.5) 

plot(plot_pgls_edmass)

# Save the plot as a rectangle with specified width and height
ggsave("../Plots/plot_pgls_edmass.png", plot = plot_pgls_edmass, width = 9, height = 7.5, units = "in")

PGLS of log eye diameter vs. log snout-vent length

# PGLS eye diameter vs. SVL
pgls_edsvl <- pgls(log10(eye_av) ~ log10(svl_av), 
               data = sal.comp, 
               lambda = "ML", #uses Maximum Likelihood estimate of lambda
               param.CI = 0.95)

# Diagnostic plots
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_edsvl)

par(mfrow = c(1, 1))

# Likelihood plot for Pagel's lambda. Solid red line indicates estimate for lambda and broken red lines indicate the 95% confidence interval
lambda.edsvl <- pgls.profile(pgls_edsvl, "lambda")
plot(lambda.edsvl)

# Print model output 
summary(pgls_edsvl)
## 
## Call:
## pgls(formula = log10(eye_av) ~ log10(svl_av), data = sal.comp, 
##     lambda = "ML", param.CI = 0.95)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0206434 -0.0040166  0.0003049  0.0037003  0.0175331 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.779
##    lower bound : 0.000, p = 5.7065e-14
##    upper bound : 1.000, p = 2.6052e-08
##    95.0% CI   : (0.542, 0.918)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)   -0.887868   0.082795 -10.724 < 2.2e-16 ***
## log10(svl_av)  0.809726   0.041312  19.600 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.005968 on 137 degrees of freedom
##   (16 observations deleted due to missingness)
## Multiple R-squared: 0.7371,  Adjusted R-squared: 0.7352 
## F-statistic: 384.2 on 1 and 137 DF,  p-value: < 2.2e-16
# Quick plot of pgls
plot(log10(eye_av) ~ log10(svl_av), data = av.sal)
abline(pgls_edsvl)

Plot of eye diameter versus snout-vent length

# Plot average eye diameter (mm) against snout-vent length (mm)
plot_pgls_edsvl <- ggplot(av.sal, aes(y=eye_av, x=svl_av, text=Genus_Species)) +
  geom_point(alpha = 0.9, aes(color = Adult_Habitat), size = 5, na.rm = TRUE) +  
  scale_color_manual(values = col_hab, name = "Adult habitat") +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(size = 26),  
    axis.text.y = element_text(size = 26),  
    axis.title.x = element_text(size = 30, color = "black", margin = margin(t = 10)),  
    axis.title.y = element_text(size = 30, color = "black", margin = margin(r = 18)),  
    axis.ticks.length = unit(0.3, "cm"), 
    panel.border = element_rect(size = 2) 
  ) +
  scale_y_log10(name = "Eye diameter (mm)", breaks = c(2, 3, 7)) +
  scale_x_log10(name = "Snout-vent length (mm)", breaks = c(30, 70, 150), limits = c(25,160)) +
  geom_abline(slope = coef(pgls_edsvl)[[2]], intercept = coef(pgls_edsvl)[[1]], linetype = "solid", size = 1.5)

plot(plot_pgls_edsvl)

# Save the plot as a rectangle with specified width and height
ggsave("../Plots/plot_pgls_edsvl.png", plot = plot_pgls_edsvl, width = 9, height = 7.5, units = "in")

PGLS of log eye diameter vs. log cornea diameter

# PGLS eye diameter vs. cornea diameter
pgls_edcd <- pgls(log10(cor_av) ~ log10(eye_av), 
               data = sal.comp, 
               lambda = "ML", #uses Maximum Likelihood estimate of lambda
               param.CI = 0.95)

# Diagnostic plots
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_edcd)

par(mfrow = c(1, 1))

# Likelihood plot for Pagel's lambda. Solid red line indicates estimate for lambda and broken red lines indicate the 95% confidence interval
lambda.edcd <- pgls.profile(pgls_edcd, "lambda")
plot(lambda.edcd)

# Print model output 
summary(pgls_edcd)
## 
## Call:
## pgls(formula = log10(cor_av) ~ log10(eye_av), data = sal.comp, 
##     lambda = "ML", param.CI = 0.95)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0078114 -0.0015818 -0.0000318  0.0016618  0.0076100 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.514
##    lower bound : 0.000, p = 0.00035263
##    upper bound : 1.000, p = 8.8929e-14
##    95.0% CI   : (0.132, 0.817)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)   -0.141399   0.017157 -8.2416 1.223e-13 ***
## log10(eye_av)  0.980989   0.021445 45.7444 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.002427 on 137 degrees of freedom
##   (16 observations deleted due to missingness)
## Multiple R-squared: 0.9386,  Adjusted R-squared: 0.9381 
## F-statistic:  2093 on 1 and 137 DF,  p-value: < 2.2e-16
# Quick plot of pgls
plot(log10(cor_av) ~ log10(eye_av), data = av.sal)
abline(pgls_edcd)

Plot of eye diameter versus cornea diamter

# Plot average eye diameter (mm) against cornea diameter (mm)
plot_pgls_edcd <- ggplot(av.sal, aes(y=cor_av, x=eye_av, text=Genus_Species)) +
  geom_point(alpha = 0.9, aes(color = Adult_Habitat), size = 5, na.rm = TRUE) +  
  scale_color_manual(values = col_hab, name = "Adult habitat") +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(size = 26), 
    axis.text.y = element_text(size = 26), 
    axis.title.x = element_text(size = 30, color = "black", margin = margin(t = 10)), 
    axis.title.y = element_text(size = 30, color = "black", margin = margin(r = 18)), 
    axis.ticks.length = unit(0.3, "cm"),  
    panel.border = element_rect(size = 2) 
  ) +
  scale_y_log10(name = "Cornea diameter (mm)") +
  scale_x_log10(name = "Eye diameter (mm)", breaks = c(2, 3, 7)) +
  geom_abline(slope = coef(pgls_edcd)[[2]], intercept = coef(pgls_edcd)[[1]], linetype = "solid", size = 1.5)

plot(plot_pgls_edcd)

# Save the plot as a rectangle with specified width and height
ggsave("../Plots/plot_pgls_edcd.png", plot = plot_pgls_edcd, width = 9, height = 7.5, units = "in")

PGLS of log cornea diameter vs. log cube root of mass

# PGLS cornea diameter vs. cube root of mass
pgls_cdmass <- pgls(log10(cor_av) ~ log10(rootmass_av),
               data = sal.comp, 
               lambda = "ML", #uses Maximum Likelihood estimate of lambda
               param.CI = 0.95)

# Diagnostic plots
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_cdmass)

par(mfrow = c(1, 1))

# Likelihood plot for Pagel's lambda. Solid red line indicates estimate for lambda and broken red lines indcaite the 95% confidence interval
lambda.cdmass <- pgls.profile(pgls_cdmass, "lambda")
plot(lambda.cdmass)

# Print model output 
summary(pgls_cdmass)
## 
## Call:
## pgls(formula = log10(cor_av) ~ log10(rootmass_av), data = sal.comp, 
##     lambda = "ML", param.CI = 0.95)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0298507 -0.0045290  0.0006747  0.0073960  0.0295073 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.953
##    lower bound : 0.000, p = < 2.22e-16
##    upper bound : 1.000, p = 0.00068248
##    95.0% CI   : (0.885, 0.988)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        0.154229   0.059419  2.5956  0.01036 *  
## log10(rootmass_av) 0.654704   0.042908 15.2584  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01015 on 153 degrees of freedom
## Multiple R-squared: 0.6034,  Adjusted R-squared: 0.6008 
## F-statistic: 232.8 on 1 and 153 DF,  p-value: < 2.2e-16
# Quick plot of pgls
plot(log10(cor_av) ~ log10(rootmass_av), data = av.sal)
abline(pgls_cdmass)

Plot of cornea diameter versus cube root of mass

# Plot average cornea diameter (mm) against cube root of mass (g)
plot_pgls_cdmass <- ggplot(av.sal, aes(y=cor_av, x=rootmass_av, text=Genus_Species)) +
  geom_point(alpha = 0.9, aes(color = Adult_Habitat), size = 5, na.rm = TRUE) +  
  scale_color_manual(values = col_hab, name = "Adult habitat") +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(size = 26),  
    axis.text.y = element_text(size = 26), 
    axis.title.x = element_text(size = 30, color = "black", margin = margin(t = 10)), 
    axis.title.y = element_text(size = 30, color = "black", margin = margin(r = 18)), 
    axis.ticks.length = unit(0.3, "cm"),  
    panel.border = element_rect(size = 2) 
  ) +
  scale_y_log10(name = "Cornea diameter (mm)") +
  scale_x_log10(name = "Cube root of mass (g)", breaks = c(1, 3, 5), limits = c(0.5, 5.5)) +
  geom_abline(slope = coef(pgls_cdmass)[[2]], intercept = coef(pgls_cdmass)[[1]], linetype = "solid", size = 1.5)

plot(plot_pgls_cdmass)

# Save the plot as a rectangle with specified width and height
ggsave("../Plots/plot_pgls_cdmass.png", plot = plot_pgls_cdmass, width = 9, height = 7.5, units = "in")

PGLS of log cornea diameter vs. log snout-vent length

# PGLS cornea diameter vs. SVL
pgls_cdsvl <- pgls(log10(cor_av) ~ log10(svl_av),
               data = sal.comp, 
               lambda = "ML", #uses Maximum Likelihood estimate of lambda
               param.CI = 0.95)

# Diagnostic plots
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_cdsvl)

par(mfrow = c(1, 1))

# Likelihood plot for Pagel's lambda. Solid red line indicates estimate for lambda and broken red lines indcaite the 95% confidence interval
lambda.cdsvl <- pgls.profile(pgls_cdsvl, "lambda")
plot(lambda.cdsvl)

# Print model output 
summary(pgls_cdsvl)
## 
## Call:
## pgls(formula = log10(cor_av) ~ log10(svl_av), data = sal.comp, 
##     lambda = "ML", param.CI = 0.95)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.034536 -0.004406  0.000953  0.007409  0.029583 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.958
##    lower bound : 0.000, p = < 2.22e-16
##    upper bound : 1.000, p = 8.8349e-06
##    95.0% CI   : (0.904, 0.985)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)   -0.809986   0.104842 -7.7258 1.363e-12 ***
## log10(svl_av)  0.607606   0.042201 14.3980 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01066 on 153 degrees of freedom
## Multiple R-squared: 0.5754,  Adjusted R-squared: 0.5726 
## F-statistic: 207.3 on 1 and 153 DF,  p-value: < 2.2e-16
# Quick plot of pgls
plot(log10(cor_av) ~ log10(svl_av), data = av.sal)
abline(pgls_cdsvl)

Plot of cornea diameter versus snout-vent length

# Plot average cornea diameter (mm) against snout-vent length (mm)
plot_pgls_cdsvl <- ggplot(av.sal, aes(y=cor_av, x=svl_av, text=Genus_Species)) +
  geom_point(alpha = 0.9, aes(color = Adult_Habitat), size = 5, na.rm = TRUE) +  
  scale_color_manual(values = col_hab, name = "Adult habitat") +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(size = 26), 
    axis.text.y = element_text(size = 26), 
    axis.title.x = element_text(size = 30, color = "black", margin = margin(t = 10)),  
    axis.title.y = element_text(size = 30, color = "black", margin = margin(r = 18)),  
    axis.ticks.length = unit(0.3, "cm"),  
    panel.border = element_rect(size = 2) 
  ) +
  scale_y_log10(name = "Cornea diameter (mm)") +
  scale_x_log10(name = "Snout-vent length (mm)", breaks = c(30, 70, 150), limits = c(25,160)) +
  geom_abline(slope = coef(pgls_cdsvl)[[2]], intercept = coef(pgls_cdsvl)[[1]], linetype = "solid", size = 1.5)

plot(plot_pgls_cdsvl)

# Save the plot as a rectangle with specified width and height
ggsave("../Plots/plot_pgls_cdsvl.png", plot = plot_pgls_cdsvl, width = 9, height = 7.5, units = "in")

PGLS of log snout-vent length vs. log cube root of mass

# PGLS svl vs. cube root of mass
pgls_svlmass <- pgls(log10(svl_av) ~ log10(rootmass_av),
               data = sal.comp, 
               lambda = "ML", #uses Maximum Likelihood estimate of lambda
               param.CI = 0.95)

# Diagnostic plots
par(mar = c(4,4,2,2))
par(mfrow = c(2, 2))
plot(pgls_svlmass)

par(mfrow = c(1, 1))

# Likelihood plot for Pagel's lambda. Solid red line indicates estimate for lambda and broken red lines indcaite the 95% confidence interval
lambda.svlmass <- pgls.profile(pgls_svlmass, "lambda")
plot(lambda.svlmass)

# Print model output 
summary(pgls_svlmass)
## 
## Call:
## pgls(formula = log10(svl_av) ~ log10(rootmass_av), data = sal.comp, 
##     lambda = "ML", param.CI = 0.95)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0173001 -0.0037582 -0.0003353  0.0026855  0.0247753 
## 
## Branch length transformations:
## 
## kappa  [Fix]  : 1.000
## lambda [ ML]  : 0.921
##    lower bound : 0.000, p = < 2.22e-16
##    upper bound : 1.000, p = 1.8264e-12
##    95.0% CI   : (0.847, 0.963)
## delta  [Fix]  : 1.000
## 
## Coefficients:
##                    Estimate Std. Error t value  Pr(>|t|)    
## (Intercept)        1.620475   0.033299  48.664 < 2.2e-16 ***
## log10(rootmass_av) 0.996038   0.026331  37.827 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.005711 on 153 degrees of freedom
## Multiple R-squared: 0.9034,  Adjusted R-squared: 0.9028 
## F-statistic:  1431 on 1 and 153 DF,  p-value: < 2.2e-16
# Quick plot of pgls
plot(log10(svl_av) ~ log10(rootmass_av), data = av.sal)
abline(pgls_svlmass)

Plot of snout-vent length versus cube root of mass

# Plot average snout-vent length (mm) against cube root of mass (g)
plot_pgls_svlmass <- ggplot(av.sal, aes(y=svl_av, x=rootmass_av, text=Genus_Species)) +
  geom_point(alpha = 0.9, aes(color = Adult_Habitat), size = 5, na.rm = TRUE) +  
  scale_color_manual(values = col_hab, name = "Adult habitat") +
  theme_bw() +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    axis.text.x = element_text(size = 26.3), 
    axis.text.y = element_text(size = 26.3), 
    axis.title.x = element_text(size = 30.3, color = "black", margin = margin(t = 10)),  
    axis.title.y = element_text(size = 30.3, color = "black", margin = margin(r = 12)),  
    axis.ticks.length = unit(0.3, "cm"),  
    panel.border = element_rect(size = 2)  
  ) +
  scale_y_log10(name = "Snout-vent length (mm)", breaks = c(30, 70, 150), limits = c(25,160)) +
  scale_x_log10(name = "Cube root of mass (g)", breaks = c(1, 3, 5), limits = c(0.5,5.5)) +
  geom_abline(slope = coef(pgls_svlmass)[[2]], intercept = coef(pgls_svlmass)[[1]], linetype = "solid", size = 1.5)

plot(plot_pgls_svlmass)

# Save the plot as a rectangle with specified width and height
ggsave("../Plots/plot_pgls_svlmass.png", plot = plot_pgls_svlmass, width = 9.4, height = 7.5, units = "in")